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' Abstract. I provide a pedagogic review of adaptive mesh refinement (AMR) radiation hydro- 
dynamics (RHD) methods and codes used in simulations of star formation, at a level suitable 
for researchers who are not computational experts. I begin with a brief overview of the types of 
£jq' RHD processes that are most important to star formation, and then I formally introduce the 
equations of RHD and the approximations one uses to render them computationally tractable. 
I discuss strategies for solving these approximate equations on adaptive grids, with particular 
emphasis on identifying the main advantages and disadvantages of various approximations and 

' numerical approaches. Finally, I conclude by discussing areas ripe for improvement. 
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_J ■ 1. Introduction 

^ ■ 

O |. While gravity is the dominant force in star formation, radiative processes are crucial 
as well. Most basically, radiation removes energy, allowing a collapsing cloud to maintain 
5^ ■ a nearly constant, low temperature as its density and gravitational binding energy rise 
by many orders of magnitude. Radiative cooling is what makes star formation a dynamic 
rather than a quasi-static process, a point understood quite early (see the review by 
lHayashil 1966f ). At densities below ~ 10 4 ~ 5 H atoms cm -3 , cooling occurs primarily via 



collisionally-excited atomic or molecular line emission, while at higher densities dust 
and gas become thermally well-coupled v ia collisions, and thermal radiation by dust 
grains is the dominant cooling process (e.g. IGoldsmith 2001 ). At yet higher densities the 
coupled dust-gas fluid becomes optically thick t o its own coo ling radiation, causing the 



temperature to rise in a quasi-adiabatic fashion ( Larson! 19691 ) 



While one-dimensional simulations of this collapse can include quite sophisti cated 



1998), the 



treatments of radiative transfer between gas parcels (e.g. iMasunaga et al. 
computational difficulty of handling RHD in three dimensions led to a long period in 

t-H . which most 3D simulations of star formation simply dispensed with radiation entirely, 
choosing instead to represent the cooling processes either by a cooling function that re- 

• »"H . moves energy at a rate dependent only on local gas properties, or by an even simpler 
equation of state that prescribes the gas temperature as a function of density. However, 
this approach has a major flaw: it does not allow energy exchange between parcels of 
gas, or between stars and diffuse gas. The last few years have seen a number of studies 
pointing out that this energy exchange is crucial to the dynamics of star formation, and 
that it can be dealt with only by RHD simulations. 

The most important radiative transfer processes can be broken into three rough cat- 
egories: thermal feedback, in which collapsing gas and stars heat the gas and thereby 
change its pressure; force feedback, in which radiation exerts forces on the gas that alter 
its motion; and chemical feedback, in which radiation changes the chemical state of the 
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gas (e.g. by ionizing it), and this chemical change somehow affects the dynamics. Each 
of these processes has been the obje ct of intense s tudy in the last few years. 

With regard to thermal feedback, lLarson pointed out that how gas fragments, 

and thus the stellar mass distribution that resu lts from fragment atio n, is extremely sensi- 
tive to how the temperature varies with density. iKrumholzl ( 20061 ) and Krumholz fc McKee 
(200$]) showed that one consequence of this sensitivity is that radiation produced by the 
first few stars to form in a given region can heat the gas around them, reducing the ability 
of that gas to fragment, favoring monolithic collapse to massiv e stars and suppressing 
formation of low mass objects. Subsequent numerical si mulations (jKrumholz et al 1 l2007aL 
2010t lBatell2009t lOffner et all 120091: lUrban et ai]|2010l) using a variety of methods have 



confirmed this effect. 

Force feedback becomes important in the context of massive stars, which produce such 
high luminosities that the radiative force they exert on dusty interstellar gas can exceed 
their gravitational force. A number of authors have argue d based on one-dimensional 
models that this process sets an u pper limit on stellar masses jlKahnlll974HYorke fc Kruegell 
1977t IWolfire fe Cassineililll987h. More recent work in two (INakanolll989l: iNakano et al.l 
19951 Lfiiina fc Adamslll996l lYorke fc Bodenheimejll999t lYorke fe Sonnhalterll2002l also 



see the con tribution by Kuiper et al. in these proceedings) and three ([Krumholz et al 



2007all2009h dimensions instead suggests that force feedback does not limit stellar masses. 



Deciding the question requires radiation-hydrodynamic simulations. 

The most important type of chemical feedback is photoionization, which raises the gas 
from ~10Kto~10 4 K. This causes a nu mber of important dynamic effects, includ- 
ing ejecting mass from star-for ming clouds (Dale et al.| 2005 : Peters et al. 2010l ). alter- 
ing the magnetic field str ucture (jKrumholz et al.ll2007d l. and driving turbulent motions 
( Gritschneder et al. 2009). Radiation also drives other chemical processes, but these are 
usually not important for dynamics, and thus may be handled by post-processing simula- 
tions. In this review I limit myself to radiative effects that are dynamically important and 
must therefore be simulated in tandem with the hydrodynamic evolution of the system. 



2. The Equations of Radiation Hydrodynamics 

Th e equations of radiation hydrodynamics (RHD) in conservation form read ([Mihalas fc Weibel-Mihalas 
19991) 




pv:v + P | = ( G | , (2.1) 
(pe + P)v 

where /?, v, e, and P are the gas density, velocity, specific energy, and pressure. The 
source term on the right hand side represents the rate at which radiation transfers mo- 
mentum and energy to the gas. For simplicity I have omitted terms describing gravity 
and magnetic fields, which would appear as additional sources on the right hand side. 

The rate at which radiation transfers energy and momentum to the gas is described 
by the radiation four-force vector (G°, G), 



G° 
G 



-- J dvj dSl [/s„(n)Z tf -t) v (n)] ( M (2.2) 

where n is a unit vector, I u is the radiation intensity at frequency v in direction n, and k„ 
and r\ v are the extinction coefficient and emissivity of the g function of frequency 
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v and direction n. Finally, the intensity is governed by the transfer equation, 

~^-I v +n-VI l , = r) v (n)-K„(iL)I l ,. (2.3) 
c ot 

For simplicity I have omitted terms describing scattering, which is generally not impor- 
tant for the radiation-hydrodynamics of star formation. 

The equations of RHD arc, like the ordinary equations of fluid dynamics, characterized 
by dimensionless numbers. The two most important ones for radiation-hydrodynmics are 
(3 ~ v/c, where v is the characteristic value of v, and r ~ L/X p , where L is the size- 
scale of the system and A p is the photon mean free path. The first ratio /3 characterizes 
how relativistic the system is. In writing equation (|2.ip in non-relativistic form, we have 
already assumed that j3 <C 1. The second, r, characterizes how optically thick it is. 
Systems with t<1 are described as being in the streaming limit, and are characterized 
by weak matter-radiation interaction. Those with r»l are in the diffusion limit and 
have strong-matter radiation coupling. Both limits appears in star formation problems. 



3. Approximations and Solution Methods 

The full system of equations (|2.1[) and (|2.3p is seven-dimensional, with quantities vary- 
ing in time, space (3 dimensions), direction (2 dimensions), and frequency. Unfortunately, 
this makes full numerical solution prohibitively expensive even on modern supercom- 
puters. We are therefore reduced to approximations. There are two broad classes of 
approximation in common use in star formation simulations: moment methods and char- 
acteristic methods. A third category, Monte Carlo methods, has been used extensively 
in post-processing radiative transfer calculations, but has not been used extensively in 
radiation-hydrodynamic simulations. I will not discuss it in detail. 

3.1. Codes 

In the discussion that follows I introduce the most common methods used in AMR 
RHD methods, and in Table [T] I summarize which method(s) are used in each code. 
The codes I include in the table are as follows: Orion is the oldest and probably best- 
tested AMR RHD code, but is not publicly available, and currently lacks a magnetohy- 
drodyn amic (MHD) capability (Klein 1999t Howell &: Greenough 20031 : iKrumholz et al 
2007bh . Ramses is an AMR MHD code that was recently upgraded with a radiation 



solve r although the latter is not (as of this writing) publicly available (jFromang et al 



2006t ICommercon et alJlioirl . Flash is an open source AMR MHD code with exten- 
sive c hemistry capabilities and several add-on modules that handle r adiation in different 
ways ( Frvxell et al. l2000t Riikhorst et al. 20061: Peters et al. 2010l) . Enzo is an open 
source AMR HD code (an MHD version exi sts but is not public) with cosmology capa- 



bilities and two different radiation methods (Abel fc Wandelt 2002; Norman et al. 2008 



Reynolds et al.1 2009). P luto is an open source AMR MHD code with a newly-developed 
radiation solver (|Mignone et al. 2007 : Kuiper et al. 2010l ). There are also a number of 
fixed, nested grid codes in use, but I will not mention these by name. 

3.2. Moment Methods 

3.2.1. Basic Theory of Moment Methods 

The basic idea of a moment method is to take moments of the transfer equation, in 
exact analogy to the Chapman-Enskog procedure used to derive the equations of fluid 
dynamics from the kinetic theory of gases. To take the zeroth moment we integrate the 
transfer equation over all directions; to obtain the first moment we multiply both sides 
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Moment methods Characteristic methods 



Code 


v resolution 


Frame 


Code 


v resolution 


Ray scheme 


Orion 
Ramses 
Enzo 
Pluto 


MG, 2T 
2T 
2T 
IT 


Mixed 
Comoving 
Comoving 
Comoving 


Orion 
Flash 
Enzo 
Pluto 


IF 

MG 
IF 
MG 


HEALpix 
Hybrid characteristics 
HEALpix 
Sphere 





Table 1. Summary of Codes. See text for explanation of fields. All moment codes use first 
order closure (flux- limited diffusion). Enzo and Orion appear in both lists because they can 
operate in characteristic or moment mode. Pluto appears ion both lists because it uses a hybrid 
characteristic- moment method (see text for details). Orion has both a MG and a 2T moment 
method. 

by the unit vector n and then integrate; for the second moment we multiply by the rank 
two tensor n : n and integrate, and so forth. As in the analogous fluid case, the procedure 
yields an exact solution if carried out to infinitely many orders, but one instead makes an 
approximation by truncating the procedure after finitely many orders. Usually "finitely 
many" here means one or two, and the first two moments of equation (|2.3j) are 

d ( E \ „ ( F \ / cG° 



at V f/c 2 j T v " V n j y G 

where 



(3.1) 



E=- I dv I dm u (3.2) 



c 



dv J dttl v n (3.3) 

n=iy dv J dni u ii:n (3.4) 

are the first three moments of the radiation intensity: the radiation energy density, ra- 
diation flux, and radiation pressure tensor. In these equations the direction-dependence 
in the transfer equation is removed and the dimensionality of the problem is lowered by 
two, but at the price of introducing the radiation pressure tensor, for which we do not 
have an equation (since we have not expanded the transfer equation to the next order) 
and must instead make an approximation. 

Because moment methods necessarily involve smoothing the angular dependence of the 
radiation field, they are best suited to representing diffuse, smooth radiation fields. This 
makes them ideal for handling thermal and force feedback, which tend to be dominated 
by diffuse infrared light re-radiated by dust grains. They are less suited for chemical 
feedback, which is usually dominated by direct, highly beamed radiation from a handful 
of stellar sources. A further advantage of moment methods is that their computational 
cost is independent of the number of sources, and usually scales only as N or NlogN, 
where N is the number of cells. They also parallelize fairly easily. 

3.2.2. Types of Moment Methods 

Moment methods can be classified based on several design choices that are made in 
their construction, all of which involve some tradeoff between computational expense and 
accuracy. See Table Q] for a list of the design choices made in each of the popular codes. 

Closure order. The first choice is whether to retain both of the first two moment 
equations, or to retain only the zeroth moment and adopt a closure approximation for 
the radiation flux; all the moment methods currently used in practical star formation 
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simulations make the latter choice, using a closure known as flux-limited diffusion (FLD). 
In a medium where the photon mean free path is small compared to the size of the system, 
Eddington showed that the radiation flux in the fluid rest frame approaches 

F = ~£-VE, (3.5) 

where kr = J dv (dB„ / dT) / J dv K~ 1 (dB u /dT) is the Rosseland mean opacity, B„{T) 
is the Planck function. This is known as the diffusion approximation. While it is quite 
good in optically thick media, it breaks down in optically thin regions because as kr — > 
the signal speed approaches infinity rather than properly limiting to c. The FLD approx- 
imation is an attempt to fix this problem by instead setting the radiation flux to 

F = -— V£ , (3.6) 

where A is a function that has the property that A — > 1/3 in optically thick regions and 
A — > KR_E /Vi?o in optically thin regions, so that |F | — > cE Q , its correct limiting value. 
Many choices for the function A are possible; the most common one in astrophysical 



applications is the iLevermore fc Pomraning (1981) limiter, A(i?) = R 1 (cothi? — R 1 ), 



where R = Vi?o|/ (%£o)- However, all of these limiters are of unknown accuracy in the 
intermediate optical depth regime, and all FLD methods suffer from the problem that 
they discard information about the directionality and momentum content of the radiation 
field. Higher order moment methods exist and can solve some of these problems (e.g. 



variable tensor Eddington fac tor methods, lHaves fc Normanll2003l and the Ml closure 



method, iGonzalez et al.l I2007T) . but none have yet proven cheap and robust enough for 
use in practical AMR star formation simulations. 

Frequency resolution. The second choice in setting up an RHD moment method is the 
level of resolution in frequency. Most accurate is the multigroup (MG) method, in which 
one discretizes equation (|3.1|) in frequency, solving one version of the equation for each 
frequency bin. All the bins are coupled via the matter temperature, which affects the 
extinction and emissivity and thus the radiation four- force vector (equation I2.2p . Next 
most accurate is the 2T method, in which one takes the radiation spectrum to be a Planck 
function characterized by a single temperature T ra d (thus removing a dimension from the 
problem), but allows T rac j to be different than the gas temperature T gas . Simplest of all 
is the IT method, in which one assumes T ra( j = Tg as , allowing a further simplification of 
the equations. As one might expect, this choice involves a tradeoff between accuracy and 
cost; the IT method is cheapest but badly misestimates both the radiation spectrum 
and T gas in the streaming regime. In comparison 2T methods still produce incorrect 
radiation spectra in the streaming regime, but make much smaller errors in the matter 
temperature. They are intermediate in cost. MG provides a good representation of both 
the spectrum and the matter temperature but is most expensive. 

Choice of frame. The third design choice is whether to formulate the equations in 
the comoving frame or using mixed frames. Extinction k and emissivity rj are simple, 
isotropic functions in the frame comoving with the fluid, but in other frames they con- 
tain complex directional dependence as a result of relativistic boosting. One might think 
that this effect is small in non-relativist ic flows, but ignoring it is equivalent to neglecting 
the work done by radiation on the gas (|Mihalas fc Kleinlll982h . which is obviously unac- 
ceptable if radiation forces are non- negligible. To avoid complex velocity dependence in 
rj and k, it is highly desirable to write them in the comoving frame. The simplest choice 
aft er that is to write the r adiation quantities (E, F) in the comoving frame as well (e.g. 
see Mihalas &: Kleinlll982 ). but this carries a significant cost. Since the comoving frame 
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is non-inertial in any system with non-constant fluid velocity, the comoving radiation 
energy is not a conserved quantity, and thus comoving-frame codes cannot be exactly 
conservative. Moreover, every time the resolution changes and the velocity field is refined, 
the reference frame change as well, introducing errors in conservation that are likely to 
accumulate with increasing refinement. The alternative is a mixed- frame formulation in 
which one writes the radiation quantities in the inertial lab frame. Deriving this formula- 
tion is tricky, since one must carefully account for the Lorentz transformations between 
frames, but the resulting equations are explicitly cons ervative, and can be dis cretized in 



a manner that conserves energy to machine precision (jKrumholz et al.ll2007bl ) 



3.3. Characteristic Methods 

In characteristic methods one solves the transfer equation (|2.3I) directly, but only along 
selected rays. Given this solution, one can compute the radiation four-force vector (equa- 
tion I2.2[) directly. This provides much greater accuracy in computing the radiation in- 
tensity along those rays, but at the cost of neglecting all other rays. For this reason it is 
best used for chemical feedback, which tends to be dominated by the rays coming from 
a small number of sources. It is not well-suited to handling diffuse radiation fields, and 
simulations based on this technique are generally no better than non-radiative simula- 
tions when it comes to handling effects like fragmentation being altered by a diffuse IR 
radiation field. The cost of these methods scales as the number of cells times the number 
of sources, but with a coefficient that is generally smaller than for moment methods. 
Thus these methods are cheaper than moment methods when the number of sources and 
rays is small, but become impractically expensive if the number of sources is even a small 
fraction of the number of cells. As with moment methods, there are several design choices 
to be made in setting up a characteristic method. 

Frequency resolution. As in the moment method case, one can choose either to retain 
the frequency-dependence in the transfer equation by using multiple frequency groups 
(the MG method), or one can integrate over frequency. In this case one does not generally 
assign a temperature to the radiation; instead, since for chemical feedback one cares about 
photons only above a certain energy, the usual approximation is to adopt an average 
photon energy. This is the single- frequency (IF) approximation. 

Ray-drawing scheme. There are also a number of schemes for drawing rays from the 
point source(s). The simplest is the spherical ray technique: one uses spherical coordi- 
nates, requires that the source be at the origin, and simply draws rays aligned with the 
computational grid. This parallelizes almost perfectly and is very simple to code, but 
is obviously unsuited to any problem with multiple stellar sources or where the loca- 
tions at which stars form is not known a priori. The ot her two schemes in wide use are 



HEALPix-based adap tive trees (jAbel fc Wandeltl 120021) and the hybrid characteristics 



(|Riikhorst et al.ll2006l ). both of which operate on top of an underlying Cartesian grid. 
The HEALPix scheme is based on the Hierarchical Equal Area isoLatitude Pixelization 
of the sphere, which divides the sphere into a equal area pixels that can be subdivided 
indefinitely, allowing the angular resolution to adapt to match that of the underlying 
grid. The hybrid characteristics scheme works by combining a method of short char- 
acteristics within individual AMR grids with a method of long characteristics between 
grids. Neither suffer from the limitations of the spherical ray method. 

3.4. Hybrid Methods 

Hybrid methods attempt to combine the best features of characteristic and moment 
methods. Characteristics work well for the highly beamed radiation coming directly from 
a star or some other point source. However, once this radiation is absorbed, following its 
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re-emission with characteristics becomes unreasonably expensive, and a moment method 
is a far better choice. Th e underlying idea of a hybrid method, first used in multidi- 
mensi onal simulations b y IMurrav et al. and recently implemented in the Pluto 

code ( Kuiper et al. 2010h , is to use characteristics for the "first absorption" , then switch 



to a moment method. One does this by performing a characteristic trace from the star 
or stars to determine the rate of energy and momentum input into each cell by direct 
radiation. One adds this as a source term on the right hand side of the moment equation 
(|3.1|) . then solves as one would in a pure moment method. The computationally-cheap 
characteristic step can use MG, while the more expensive moment solve uses the IT or 
2T approximation. Since the first absorption is generally the part of the problem where 
the frequency dependence is most important, this method achieves some of the accuracy 
of a fully frequency-dependent calculation at significantly lower cost. 



4. Future Directions 

I close this review with a brief discussion of possible future directions for AMR RHD 
in the star formation context. One such likely direction comes from combining realistic 
treatments of thermal, force, and chemical feedback into a single simulation, rather than 
treating only one or two of them at a time as in current simulations. Since the first two 
effects are most easily handled by moment methods and the third by characteristic meth- 
ods, hybrid techniques are the natural solution. To be competitive in handling thermal 
and force feedback, however, these methods will have to be linked to more advanced 
moment methods than has been attempted before. Nonetheless, a hybrid method com- 
bining frequency-dependent characteristic tracing with a 2T, mixed frame FLD method 
is a straightforward extension of current techniques, and seems a logical place to start. 

Further in the future, second order moment methods are likely to become important. 
Developing them to the point where they can run reliably in a parallel environment on 
adaptive grids will be a major algorithmic undertaking, one that is likely to require the 
assistance of computer scientists and applied mathematicians in addition to astronomers. 
If successful, these methods will be able to handle both beamed and diffuse radiation 
fields within a single framework, and will remove many of the limitations of the FLD 
approximation. Monte Carlo methods are another candidate to replace the current dom- 
inance of FLD and characteristic methods. They are very computationally costly, since 
many photons are required to suppress Poisson noise, but they have the advantage of near 
perfect parallelization. Since the future of supercomputing seems to be heading toward 
the development of ever-larger numbers of processors, each of which is no faster than the 
processors of the previous generation, this may prove to be a decisive advantage. 

Finally, a caveat: for all the limitations of our current RHD methods, it is not entirely 
obvious that our numerical techniques are the limiting factor in the accuracy of our sim- 
ulations. The main agent for coupling radiation and gas at the high densities where stars 
form is dust. Our knowledge of dust grain properties, such as sublimation temperatures 
and grain size distributions in dense environments, is quite limited, and appears likely to 
advance less quickly than either computer power or algorithmic technique. It may be the 
case in the future that our knowledge of dust becomes the limiting factor in the accuracy 
of RHD simulations. 
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